Maximizing as minimum
Adapted from: (Floudas et al., 1999; Section 4.10), (Laurent, 2008; Example 6.23) and (Lasserre, 2009; Table 5.1)
Introduction
Consider the polynomial optimization problem from (Floudas et al., 1999; Section 4.10) of minimizing the linear function $-x_1 - x_2$ over the basic semialgebraic set defined by the inequalities $x_2 \le 2x_1^4 - 8x_1^3 + 8x_1^2 + 2$, $x_2 \le 4x_1^4 - 32x_1^3 + 88x_1^2 - 96x_1 + 36$ and the box constraints $0 \le x_1 \le 3$ and $0 \le x_2 \le 4$,
using DynamicPolynomials
@polyvar x[1:2]
p = -sum(x)
using SumOfSquares
f1 = 2x[1]^4 - 8x[1]^3 + 8x[1]^2 + 2
f2 = 4x[1]^4 - 32x[1]^3 + 88x[1]^2 - 96x[1] + 36
K = @set x[1] >= 0 && x[1] <= 3 && x[2] >= 0 && x[2] <= 4 && x[2] <= f1 && x[2] <= f2Basic semialgebraic Set defined by no equality
6 inequalities
x[1] ≥ 0
3 - x[1] ≥ 0
x[2] ≥ 0
4 - x[2] ≥ 0
2 - x[2] + 8*x[1]^2 - 8*x[1]^3 + 2*x[1]^4 ≥ 0
36 - x[2] - 96*x[1] + 88*x[1]^2 - 32*x[1]^3 + 4*x[1]^4 ≥ 0
As we can observe below, the bounds on x[2] could be dropped and optimization problem is equivalent to the maximization of min(f1, f2) between 0 and 3.
xs = range(0, stop = 3, length = 100)
using Plots
plot(xs, f1.(xs), label = "f1")
plot!(xs, f2.(xs), label = "f2")
plot!(xs, 4 * ones(length(xs)), label = nothing)
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.OptimizerClarabel.MOIwrapper.OptimizerA Sum-of-Squares certificate that $p \ge \alpha$ over the domain S, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
endsolve (generic function with 1 method)The first level of the hierarchy gives a lower bound of -7`
model4 = solve(4)-------------------------------------------------------------
Clarabel.jl v0.7.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 31
constraints = 40
nnz(P) = 0
nnz(A) = 73
cones (total) = 6
: Zero = 1, numel = 10
: PSDTriangle = 5, numel = (6,6,6,6,6)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 1.2590e-01 1.2590e-01 8.33e-17 6.99e-01 6.62e-01 1.00e+00 1.51e+00 ------
1 4.8332e-01 5.0333e-01 2.00e-02 1.40e-01 9.41e-02 1.86e-01 3.32e-01 8.02e-01
2 1.3995e+00 1.5138e+00 8.17e-02 8.66e-02 2.59e-02 2.10e-01 1.08e-01 7.85e-01
3 2.7854e+00 2.9400e+00 5.55e-02 2.03e-02 3.36e-03 1.83e-01 1.58e-02 8.78e-01
4 4.0242e+00 4.2592e+00 5.84e-02 1.43e-02 1.34e-03 2.61e-01 7.20e-03 7.46e-01
5 5.4090e+00 5.5096e+00 1.86e-02 4.16e-03 3.04e-04 1.09e-01 1.81e-03 8.27e-01
6 6.2811e+00 6.3385e+00 9.13e-03 1.73e-03 1.07e-04 6.15e-02 6.79e-04 7.22e-01
7 6.9685e+00 6.9720e+00 5.12e-04 6.60e-05 3.84e-06 3.75e-03 2.52e-05 9.69e-01
8 6.9993e+00 6.9994e+00 1.29e-05 1.66e-06 9.63e-08 9.46e-05 6.33e-07 9.75e-01
9 7.0000e+00 7.0000e+00 3.48e-07 4.54e-08 2.62e-09 2.56e-06 1.72e-08 9.73e-01
10 7.0000e+00 7.0000e+00 1.18e-08 1.56e-09 8.90e-11 8.65e-08 5.85e-10 9.66e-01
11 7.0000e+00 7.0000e+00 2.34e-10 3.13e-11 1.84e-12 1.72e-09 1.17e-11 9.80e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 3.34ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -7.00000e+00
Dual objective value : -7.00000e+00
* Work counters
Solve time (sec) : 3.33579e-03
Barrier iterations : 11The second level improves the lower bound
model5 = solve(5)-------------------------------------------------------------
Clarabel.jl v0.7.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 108
constraints = 128
nnz(P) = 0
nnz(A) = 266
cones (total) = 7
: Zero = 1, numel = 21
: Nonnegative = 1, numel = 2
: PSDTriangle = 5, numel = (21,21,21,21,21)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 2.8313e-01 2.8313e-01 0.00e+00 7.53e-01 6.03e-01 1.00e+00 1.82e+00 ------
1 5.4905e-01 5.6026e-01 1.12e-02 1.28e-01 7.06e-02 1.86e-01 4.16e-01 8.05e-01
2 1.0075e+00 1.0482e+00 4.04e-02 3.51e-02 1.22e-02 8.87e-02 1.08e-01 8.48e-01
3 1.7352e+00 1.7929e+00 3.33e-02 8.17e-03 1.68e-03 7.10e-02 1.55e-02 9.20e-01
4 2.2703e+00 2.3425e+00 3.18e-02 3.66e-03 5.84e-04 7.97e-02 4.20e-03 8.20e-01
5 2.5204e+00 2.5826e+00 2.47e-02 2.51e-03 3.63e-04 6.81e-02 2.42e-03 6.21e-01
6 3.3688e+00 3.4392e+00 2.09e-02 6.50e-04 6.62e-05 7.25e-02 3.73e-04 8.78e-01
7 4.1613e+00 4.2336e+00 1.74e-02 3.75e-04 2.25e-05 7.36e-02 1.14e-04 8.80e-01
8 4.9013e+00 4.9482e+00 9.56e-03 1.66e-04 8.18e-06 4.76e-02 4.21e-05 8.78e-01
9 5.7265e+00 5.7548e+00 4.95e-03 6.35e-05 2.06e-06 2.86e-02 1.12e-05 8.50e-01
10 5.7843e+00 5.8135e+00 5.04e-03 6.09e-05 1.89e-06 2.95e-02 1.03e-05 2.28e-01
11 6.4423e+00 6.4536e+00 1.75e-03 1.28e-05 3.67e-07 1.13e-02 2.07e-06 8.34e-01
12 6.6492e+00 6.6500e+00 1.19e-04 7.56e-07 2.24e-08 7.94e-04 1.26e-07 9.69e-01
13 6.6662e+00 6.6663e+00 3.10e-06 1.97e-08 5.83e-10 2.08e-05 3.28e-09 9.74e-01
14 6.6667e+00 6.6667e+00 8.84e-08 5.65e-10 1.66e-11 5.93e-07 9.36e-11 9.72e-01
15 6.6667e+00 6.6667e+00 2.63e-09 1.69e-11 4.96e-13 1.76e-08 2.79e-12 9.70e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 8.35ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -6.66667e+00
Dual objective value : -6.66667e+00
* Work counters
Solve time (sec) : 8.34640e-03
Barrier iterations : 15The third level finds the optimal objective value as lower bound...
model7 = solve(7)-------------------------------------------------------------
Clarabel.jl v0.7.1 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 288
constraints = 323
nnz(P) = 0
nnz(A) = 739
cones (total) = 8
: Zero = 1, numel = 36
: PSDTriangle = 7, numel = (55,55,55,55,...,55)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 3.3234e-01 3.3234e-01 2.22e-16 7.75e-01 7.03e-01 1.00e+00 1.51e+00 ------
1 4.7370e-01 4.8104e-01 7.34e-03 1.20e-01 7.54e-02 1.77e-01 3.55e-01 8.03e-01
2 7.3313e-01 7.6805e-01 3.49e-02 4.16e-02 1.96e-02 1.12e-01 1.54e-01 7.48e-01
3 1.1439e+00 1.1684e+00 2.14e-02 7.53e-03 2.35e-03 4.09e-02 2.71e-02 8.65e-01
4 1.5684e+00 1.6151e+00 2.98e-02 3.00e-03 6.64e-04 5.81e-02 7.91e-03 8.51e-01
5 1.8613e+00 1.8904e+00 1.57e-02 1.41e-03 2.53e-04 3.50e-02 3.10e-03 6.74e-01
6 2.1699e+00 2.1913e+00 9.85e-03 5.14e-04 7.66e-05 2.40e-02 8.88e-04 8.03e-01
7 2.5021e+00 2.5220e+00 7.93e-03 1.85e-04 2.34e-05 2.11e-02 2.37e-04 8.01e-01
8 2.7896e+00 2.8141e+00 8.76e-03 1.17e-04 1.10e-05 2.54e-02 9.20e-05 8.08e-01
9 3.2619e+00 3.2812e+00 5.91e-03 3.96e-05 2.84e-06 1.96e-02 2.16e-05 7.88e-01
10 3.6591e+00 3.6805e+00 5.86e-03 2.80e-05 1.09e-06 2.17e-02 7.16e-06 9.60e-01
11 4.0143e+00 4.0278e+00 3.36e-03 1.31e-05 4.96e-07 1.36e-02 3.07e-06 7.92e-01
12 4.3776e+00 4.3897e+00 2.76e-03 7.68e-06 2.25e-07 1.22e-02 1.31e-06 6.54e-01
13 4.8717e+00 4.8798e+00 1.67e-03 2.67e-06 6.42e-08 8.16e-03 3.48e-07 8.11e-01
14 5.0064e+00 5.0151e+00 1.74e-03 2.22e-06 4.55e-08 8.76e-03 2.29e-07 5.98e-01
15 5.3203e+00 5.3231e+00 5.17e-04 6.99e-07 1.36e-08 2.76e-03 6.96e-08 9.90e-01
16 5.4810e+00 5.4813e+00 6.57e-05 1.07e-07 2.07e-09 3.62e-04 9.89e-09 8.52e-01
17 5.5042e+00 5.5043e+00 1.24e-05 1.45e-08 2.83e-10 6.85e-05 1.36e-09 9.58e-01
18 5.5078e+00 5.5078e+00 6.98e-07 8.14e-10 1.59e-11 3.86e-06 7.64e-11 9.44e-01
19 5.5080e+00 5.5080e+00 5.33e-08 5.62e-11 1.10e-12 2.94e-07 5.30e-12 9.84e-01
20 5.5080e+00 5.5080e+00 1.09e-09 1.63e-12 2.93e-14 6.04e-09 1.09e-13 9.79e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 41.1ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -5.50801e+00
Dual objective value : -5.50801e+00
* Work counters
Solve time (sec) : 4.11313e-02
Barrier iterations : 20...and proves it by exhibiting the minimizer.
ν7 = moment_matrix(model7[:c])
η = atomic_measure(ν7, 1e-3) # Returns nothing as the dual is not atomicAtomic measure on the variables x[1], x[2] with 1 atoms:
at [2.329520151218628, 3.178493153833861] with weight 0.9999999590079649We can indeed verify that the objective value at x_opt is equal to the lower bound.
x_opt = η.atoms[1].center
p(x_opt)-5.508013305052489We can see visualize the solution as follows:
scatter!([x_opt[1]], [x_opt[2]], markershape = :star, label = nothing)
This page was generated using Literate.jl.